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ABSTRACT 



In this paper, the stability of a dynamically condensing radiative gas layer is investigated by linear analysis. Our own time-dependent, 
self-similar solutions describing a dynamical condensing radiative gas layer are used as an unperturbed state. We consider pertur- 
bations that are both perpendicular and parallel to the direction of condensation. The transverse wave number of the perturbation is 
defined by k. For A: = 0, it is found that the condensing gas layer is unstable. However, the growth rate is too low to become nonlinear 
during dynamical condensation. For k t 0, in general, perturbation equations for constant wave number cannot be reduced to an 
eigenvalue problem due to the unsteady unperturbed state. Therefore, direct numerical integration of the perturbation equations is per- 
formed. For comparison, an eigenvalue problem neglecting the time evolution of the unperturbed state is also solved and both results 
agree well. The gas layer is unstable for all wave numbers, and the growth rate depends a little on wave number The behaviour of 
the perturbation is specified by kLcooi at the centre, where the cooling length, Lcooi. represents the length that a sound wave can travel 
during the cooling time. For fcLcooi s> 1, the perturbation grows isobarically. For kL^ooi *k 1, the perturbation grows because each 
part has a difi"erent collapse time without interaction. Since the growth rate is sufficiently high, it is not long before the perturbations 
become nonlinear during the dynamical condensation. Therefore, according to the linear analysis, the cooling layer is expected to 
split into fragments with various scales. 
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1. Introduction 

In the interstellar medium (ISM), it is well known that a clumpy 
low-temperature phase (cold neutral medium, or CNM) and 
a diffuse high-temperature phase (warm neutral medium, or 
WNM) can coexist in pressure equilibrium as a result of the 
balance of radiative cooling and heating due to external radia- 
tion fields and cosm ic rays (Field, Goldsmith & Habing, 1 19691: 
IWolfire et al.L [19951 120031) . These two phases are thermally sta- 
ble. On the other hand, gas is thermally unstable in the temper- 
ature range between two stable phases, that is, in the range 300 
■ K < r < 6000 K. The unstable gas spontaneously turns into the 
' mixture of stable CNM and WNM by thermal instability (TI). 
This instability is expected to play an important role in the struc- 
ture formation and the dynamics of the ISM, and especially, in 
' the molecular cloud formation and origin of turbulence. 

The basic properties of TI was investigated by Fieldl (1 19651) . 
who performed linear analysis of an uniform gas in thermally 
equilibri um. He derived a criterion for TI. Focusing on one fluid 
element, iBalbusI (11986 !) generalized the Field criterion when 
the cooling rate is not equal to the heating r ate. The effe ct of 
magnetic field on TI has been investigated by Fieldl (1 19651) and 
iHennebelle & Passo3 ( 12006 ). and other authors. 

Recently, many authors have used multi-dimensional nu- 
merical simulations to study the turbulent CNM forma- 
tion driven by TI. iKoyama & Inutsuka (2002) suggest that 
the turbulent CNM formation is induced by TI in a 
shock-compressed region. Analogous processes have been 
studied by many authors for a colliding flow of the 
WNM (Audit & Hennebelle', '2005*; Hennebelle & Audii l2007t 
[Heitsch et al., 2006; Vazquez-Semadeni et al., 2007^ and using 
two-fluid MHD simulation dlnoue & Inutsukal |2008|) . The un- 
balance between cooling and heating rates causes the shock- 



compressed gas layer to cool and to condense. During the cool- 
ing, these numerical simulations shows that the runaway cooling 
layer quickly fragments into many CNM clumps whose velocity 
dispersion is equal to a fraction of the sound speed of WNM, 
where CNM clumps and WNM are tightly interwoven. This 
complex structure is regarded as produced by TI and possibly 
by some other hydrody namical in stabilities, such as the nonlin- 
ear thin shell instability ( IVishniaa.1994.) . the Kelven-Helmholtz 
instability, and by corrugation insta bility of the phase transition 
layers between CNM and WNM dlnoue. Inutsuka & Kovamal 
I2006h . 

A fluid element that is compressed by a shock wave tends 
to be a layer rather than a sphere because it is only com- 
pressed in the direction perpendicular to the shock front. Once 
the fluid element enters the thermally unstable regime, the layer 
cools in a runaway fashion. In this paper, we focus on the frag- 
mentation of the runaway cooling layer. In previous studies, A 
detailed physical mechanism of the fragmentation of the run- 
away cooling layer remains poorly understood even in linear 
regime. The main reason is that it is difficult to select the un- 
perturbed state since the cooling layer evolves temporarily and 
spatially. Therefore, in previous works, the unperturbed states 
were limited to spatially uniform gas that cools isochorically 
(ISchwarz McCrav & Steinlll972l; iBurkert & Linl l200a) or iso- 



baricallv iKovama & InutsukaT 20C)0h 



llwasaki & Tsurib j ( 2008) (hereafter IT08) have recently 
found a family of self-similar (S-S) solutions describing the dy- 
namical condensation of a radiative gas layer where the cooling 
rate dominates the heating rate. This S-S solution assumed that 
the net cooling rate is a power-law function and that the heating 
rate is explicitly neglected. Although it is still too ideal, they are 
expected to be a good nonlinear one-dimensional model at least 
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in the phase during the transition from WNM to CNM. In this 
paper, we adopt the S-S solutions as a more realistic unsteady 
unperturbed state than those in previous works. We perform lin- 
ear analysis of the S-S solutions against fluctuations perpendic- 
ular, as well as parallel to, the direction of the condensation. By 
performing the linear analysis, we will have some useful insights 
when and how the cooling layer fragments. Since we focus on 
the above S-S unperturbed state, the nonlinear thin shell instabil- 
ity, Kelvin-Helmholtz instability, and the corrugation instability 
are beyond the scope of this paper. 

In Sect. 12] we formulate basic equations using a zooming co- 
ordinate. Perturbation equations are derived for the linear anal- 
ysis, with a brief review of the S-S solutions. In Sect. [3] we 
investigate the stability of the S-S solutions taking only those 
fluctuations into account that are parallel to the direction of the 
condensation. In Sect. H) we consider the fluctuations that are 
both perpendicular and parallel to the direction of the condensa- 
tion. In Sect. |5] we discuss the astrophysical implication of the 
Unear analysis and effects of the thermal conduction. Our study 
is summarized in Sect. |6] 



iHanawa & Matsumotol Il999l) . We introduce the similar zoom- 
ing coordinate since this transformation makes stability analysis 
easier as follows: 
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where w is a free parameter, a corresponds to the cooling length, 
and u - I - t/tc, fc is an epoch when the central density be- 
comes infinity. In the zooming coordinate, density Q, velocity 
V, pressure II, and sound speed X are given by 

Q(t, ^, .y) = pit, X, y)/pQ(t), V{t, ^, y) = v(f, x, y)/vo(t), 
n(T, ^, y) = Pit, X, y)/PQ(t), X(t, y) = c,(f, x, y)lv^(t), (6) 



respectively, where 
vo(/) = -io(f) 



(1 - w)fc 



(7) 



2. Formulation 

We consider a dynamically condensing radiative gas layer where 
the cooling rate dominates the heating rate. The following for- 
mula is adopted as the net cooling rate per unit volume and time: 



pi:(p, p) = A„p2-«p«-i oc p'^T" erg cm 



(1) 



In ISM, in the temperature range of T ^ 10^, the main coolant 
is Bre msstralung for the solar metallicity (' Sutherland & Dopital 
Il993h . The cooling rate of Bremsstralung is approximated well 
by that with a = 0.5. In the temperature range of2xlO^K$r:$ 
10^ K, the dominant coolant is m etal lines for the solar metallic- 
ity ([Sutherland & Dopitalll993h . The cooling rate of metal lines 
corresponds to that with a < 0. However, the cooling rate can- 
not be expressed by a single a. In the temperature range between 
CNM and WNM, 300 ^ T $ 6000, the dominant coolant is CII 
dPalgarno & McCravllI972l) . In this case, the cooling rate is ap- 
proximated by that with a ^ 0.6 as is shown in Sect. 15.1! 

Basic equations for a radiative gas are the continuity equa- 
tion. 



dp 

-£+V-(pv)^0, 
ot 

the equation of motion, 



Dv 1 

— + -VP = 0, 
Df p 

and the entropy equation 
D 
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{lnPp-^)^-y"Aop^-''P''-\ 



(2) 



(3) 



(4) 



where D/Dt = d/dt + v-V indicates the Lagrangian time deriva- 
tive. 

We take the jc-axis as the direction of the condensation 
driven by the cooling and y-axis as the transverse direction. 
Since the S-S solutions are time-dependent, it is difficult to 
perform linear analysi s in th e ordinary Cartesian coordinate, 
(t,x,y). iBouquet et alj (Il985h introduced a zooming coordi- 
nate where S-S solutions appear to be stationary (also see 



Po(0 



and 



1 



(1 - W)fc 



■/(l-oj) 



Ao 



r 



(8) 



(9) 



withyS = w(3 - 2a) - 1. 

In the zooming coordinate, the basic equations (|2|-(|4]i are 
rewritten as 



DlnQ. 
Dt 



DV 1 

— + _vn = wF, 

Dt Q 



and 
1 D 



2a) 



(In HQ-'') = -j6 - y"n^-"U"-\ 

y — I Dt T ~ 1 



(10) 



(11) 



(12) 



respectively, where the operators of time and spatial derivative 
are defined by 



7^ = ^ +(V + ^ef)- V, and V= — ,xo(t) — 
Dt Ot ay 



(13) 



respectively, where eg indicates the unit vector parallel to ^- 
direction. Supplements for the derivation of Eqs. (fT0l)- (fT2l) is 
presented in AppendixlAl 

We apply the zooming transformation only in the x-direction 
but not in the y-direction. This is because the gas contracts along 
X-axis but not along y-axis in the unperturbed state. In the or- 
dinary coordinate, the transverse scale of the perturbation is ex- 
pected to be constant with time. However, if the zooming trans- 
formation is also applied in the y-direction, the transverse scale 
of the perturbation decreases with time in the ordinary coordi- 
nate, although the unperturbed gas does not contract along the 
y-axis. Therefore, we apply the zooming transformation only in 
the x-direction. 
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2. 1 . Review of self-similar solutions 

In the zooming coordinate, steady state solutions correspond to 
S-S solutions that were derived in IT08. In this section, physical 
properties of the S-S solutions are reviewed briefly. 

The S-S solutions are specified by two parameters, a and ai. 
For convenience, instead of o), we can use a parameter rj, which 
is given by 



(2-a){l -(3-2a)w) 
1 - w 



(14) 



Using these parameters (a, rf), the time dependences of the cen- 
tral density, poo, and pressure. Poo, are given by 



Poo(f)-C^"* and P„o(f)-e"'\ 



(15) 



respectively, where ftp (//) = 77/(2-0') and □'^(t;) = (1 -77)/(l -ff). 

The S-S solutions include two asymptotic solutions. For 77 ~ 
0, the time dependences of the central density and pressure are 
given by 



poo(/) ~ const, and Poo(0 ^ h 



\/(l-a) 



(16) 



respectively. This time evolution indicates the isochoric mode. 
For 77 ~ 1, the time dependences of the central density and pres- 
sure are given by 



Poo(0 and Poo(0 ~ const.. 



(17) 



respectively. This time evolution corresponds to the isobaric 
mode. Our S-S solutions exist between the two limits, < 77 < 1. 
With the condition that the increasing rate of poo(0 is equal to the 
decreasing rate of Poo(0, or Q'p(77) = apirf), the critical value of 
77 is given by 



77eq = (2 - a)/(3 - 2a). 



(18) 



Therefore, the S-S solutions for < 77 < ?7eq and ?7eq < < 1 
are close to the isochoric and the isobaric modes, respectively. 
During S-S condensation, the central pressure should not in- 
crease with time. From Eq. (flST l. the range of the value of a is 
found to be Of < 1 . 



2.2. Perturbation equations 

Perturbation on the S-S solutions is considered. Perturbed vari- 
ables are defined by 

Q = Qo(^){l+5Q(T,^,y)), 
y, = VoiO + SV,{T,^,y), 
Vy = 6VyiT,^,y), 
and 

n = mm+mT,^,y)}, (i9) 

where subscript "0" indicates the unperturbed state. As the per- 
turbation, we consider the following Fourier mode with respect 
to y. 



and k indicates the wave number of the plane wave that prop- 
agates along y-direction. Substituting Eqs. (1% and ( |20] | into 
Eqs. (fT0b-(fT2ll and hnearizing, we get the following perturbation 
equations: 



DSQ. ddV, 

+ ^ = -(lnQo)'5y, - K(T)i6Vy, 



(21) 



D6V, x?.d6n , xf. 

+ ——=-(« + V,)6V, - ^(lnno)'(5n - ,5^), (22) 
Dt 70^ 7 



Di6Vy Xl 

= -u)i5Vy + K{T)—5n, 
Dt y 

and 

DOT mo. 



(23) 



■7- 



Dr Dt 

= -(2 - a)reo5Q - (a - 1 )reoOT - (In non^)'^^, (24) 
where D/Dt ^d/dT+ Vod/d^, 

eo = 7"'' (7 - l)"o"'n«-', and ^(t) = ^xo(t). (25) 

The time-dependent factors remain in the form of k(t) = 
A:xo(t) in (l2n i- (l24l l. because the transverse scale is not zoomed 
(see Eq.|5]l as mentioned above. The leads to a problem that the 
perturbed variables cannot be expanded in the Fourier mode with 
respect to t in general. 

3. Perturbation for A; = 

In this section, we consider the perturbation parallel to the con- 
densation, or for the case with A: = 0. In this case, since the 
time-dependent factor, k(t) - kxo{T), vanishes, the perturbed 
variables can be expanded in the Fourier mode with respect to 

T as 

6Q(T,^)^5Q(^)e''\ (26) 
By Eq. ( |26] |. the time evolution of the perturbations is given by 



6Q oc , where E = 



I - OJ 



(27) 



Substituting Eq. ( l26l l into the perturbation Eqs. (I2n i- (l24l l. one 
obtains the following ordinary differential equations: 



ddQi 1 ■ei 

-W = -[723^2 Z SQ = (5Q, 6V,, OT), 

^ v^o ^0 /=I 



(28) 



where the detailed expression of Aij is shown in Appendix |B] 
Equations ( l28T l are solved as a boundary- and eigenvalue prob- 
lem. 



3. 1. Boundary conditions 

We impose the boundary conditions at f = and at the critical 
point, f = fs, where Vq - Xq. The boundary conditions at f = 
are obtained by the asymptotic limit of the perturbed variables. 
From the regularity of the perturbed variables at f = 0, we find 
that perturbed variables should have the following asymptotic 
forms: 



liin 5Q(f) ^ 6D.0, lim 6V,{^ ^ -o-(5Qof , 



^QiT, f, y) = 6Q(t, f ) expiiky), where 6Q = (dQ, 6V, dU), (20) and 



cr + (a — 2)eoo 
limOT(f) - 7 ^ r-^^Qo, 

C^o 'o- + (a- l)yeoo 



(29) 



(30) 



where too = {2w -y6(y - l))/y 

The boundary conditions at the critical point, f - fs, are ob- 
tained from the Eqs. ( |28] ). At f = fs, the denominator of the 
righthand side becomes zero. To obtain a regular solution from 



4 



Kazunari Iwasaki and Tom Tsuribe: Fragmentation of a dynamically condensing radiative layer 




Fig. 1. Growth rate, E, as a function of r] for (a) a - - 1 .0 and (b) 0.5 for the case with the perturbation parallel to the condensation. 
For comparison, the increasing rate of the unperturbed central density, apirf), and the decreasing rate of the unperturbed central 
pressure, ap{rj), are shown by the dashed and dotted hnes, respectively. 



f = to ^ = oo, the numerator of the righthand side should 
vanish. Therefore, the boundary conditions are given by the fol- 
lowing three equations. 



0, / = 1,2,3. 



(31) 



The above three equations give only one independent condition. 

3.2. Numerical method 

Solutions of Eqs. (|28] ) have three integration constants. 
Therefore, if we set two constants (5Qo, o") and impose the 
boundary condition at ^ = the solution is completely fixed. 
In this section, our numerical method for solving Eqs. (l28T l is 
described. 

We can set SQ.Q - 1 without loss of generality. For a given cr, 
we integrate Eqs. ( |28] ) from^ = to the critical point, ^ = us- 
ing a fourth-order Runge-Kutta method. Eigenvalue, cr, is mod- 
ified until the perturbed variables satisfy the boundary condition 
at ^ = using the Newton-Raphson method. After that, we in- 
tegrate Eqs. ® up to ^ = 10"^. 

3.3. Results 

Figure[T]shows the dependence of growth rate, 2, on rj for (a)a - 
-1.0 and (b)0.5. From Fig. [1] it is seen that cr > for a wide range 
of a and 77. Therefore, the perturbation is unstable. However, cr is 
smaller than ap and ap. Therefore, the growth rate is too low to 
grow sufficiently during runaway condensation. This is consis- 
tent with the results of the one-dimensional numerical simulation 
shown in Sect. 15. II 



4. Perturbation with k ^ 

For k + 0, the perturbation equations contain the time-dependent 
factor, A-(f). Here we show that this factor, k - kx^if), is related 
to the ratio of the local cooling length at the centre to the wave 
length of the perturbation. The cooling timescale at f = is given 
by 



cool 



^00 



\ — (X> 



(32) 



where superscript "0" indicates the value at f = 0. Therefore, 
using r"^^j, the collapse time can be expressed as 



. r"(r- 1)^2 — O'T-rO'- 1 ^0 

'f^ ~ - ^■'00 ^™ 'cool- 



1 -W 



(33) 



Using Eqs. (|6]l and O, at f = 0, the sound speed at the centre is 
given by 

cjjo = Xoov°, where yOf^ = fl/(l - w) = xo(0)/(l - w). (34) 
Therefore, xo(0) is given by 

1 -0) 



xo(0) 



^00 



-r° t 

Cqo'c- 



Using Eq. (|33] |. Eq. ([35} can be written as 
xo(0) = MooL°„„, = a. 



(35) 



(36) 



where Moo = y-^'Hy - m'£-"Yf^'^\ and L« = 



"^cool 00 cool 



IS 



the cooling length at the centre at f = 0. Since the S-S solu- 
tions are invariant under the transformation ^ — > m^, Vq — > mVo, 
Qo m"^'°'"''Qo, and Do m^^^^^^'llo, with a parameter m, 
one can take m in such a way that Mqo is equal to 1 . Hereafter, 
Moo is set to unity. Therefore, the time-dependent factor, k{t), is 
expressed as 

k{t) = kxoir) = fcL°„„,f!^<'""^ = ^^ooi^"' = kUUr), (37) 
where Lcooi(t) = L^oof^''- 



4.1. Static approximation 

Before we present the fully time-dependent numerical calcula- 
tion in Sect. 14.21 we at first consider a special case in which the 
time evolutions of the unperturbed S-S solutions are slower than 
the growth of perturbation. In this situation, the time evolution 
of the unperturbed state is negligible during the growth of the 
perturbations. Therefore, we set x^) to be a constant in Eqs. ( 1211 1- 
(|24] |. This approximation is also valid for k <& I/xq where the 
term, k, is negligibly small in Eqs. (|2TI)-(|24|). We use the Fourier 
mode as 



SQ = (5Q, 6V„ Sn, i6Vy). 



(38) 
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Fig. 2. Approximate dispersion relations for a - 0.5 with (a) 77 = 0.1, (b) 0.75, and (c) 0.95, where k is the nondimensional wave 
number, and S the growth rate of the perturbation. The labels of number represent the branches of modes. 



The condition under which the static approximation is valid is 
given by 



dln^e 



din r* 



dln/c 



din f» 



(5 



(2 - Q-)(3 - 2g) - 77 
2(2-ff)(l-Q') ' 
2a)/ {2(2 - a)} for 77 : 



1 



(3 -2a)/ {2(1 -a)) for 77 = 



(39) 



where the definition of S is the same as that in Eq. ( |27] ). 
Substituting Eq. ( [38] l into the perturbation Eqs. (I2ni-(l24ll. we get 



ddQi 



y2 
•^0 



- y AijSQj, 6Q = (5Q, 6V„ dU, idV,), (40) 



^0 ;=1 



where the detailed expression of A,j is shown in AppendixlB] 
We impose the boundary conditions at ^ = and ^ - 
Since we are interested in the fragmentation of the cooling layer, 
only the even mode is investigated. For the even mode, the per- 
turbed variable should have the following asymptotic forms in 
^« 1: 

hm5Q(^) ^ 5Qo, 

iim(yy,(^) - (jy,o^, 

lim5y,(^) - (JVyo, 
and 

lim(5n(f) ^ OTo + OToif^ 



(41) 



Substituting Eqs. ( |4TI ) into Eqs. ( |40l ). we obtain the following 
relations: 



crQo + ,0 + 7c/(jyvo = 0, 

(cr - ljS)l5yyQ = OTo, 

r 



(42) 
(43) 



{cr + (a - 2)eoo) r<5i^o - {o" + (a - l)reoo) <Jno = 0, 



(44) 



and 



(2y8+l 



■ CO + cr)6V,-o + l-^STloi = 2j8w(l - a)(6Uo - 5Qo).(45) 

r 



The boundary condition at ^ = is derived in the same way in 
Sect. 13. II and we obtain two independent conditions. Numerical 
method for solving Eqs. ( l40l l is the same as that in Sect. 13.21 

Figures |2] shows an approximate dispersion relation for 
a = 0.5 with (a)77 = 0.95, (b) 0.75, and (c) 0.10. In Fig. |2l 
one can see several branches labelled by numbers. The most 
unstable branch is labelled by (1) for large k limit, and (2) for 
small K limit. Each branch is explained below. 




Fig. 3. Eigenfunctions for 77 = 0.75 and tc = 8. The filled cir- 
cles indicate the values at the critical point. The eigenvalue, 2, is 
equal to 1.196. 



(1) The isobaric mode 

The branch (1) corresponds to the most unstable mode for k » I, 
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(a)a=0.5 (b)a=0.8 
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Fig. 4. Growth rate as a function of rj for (a)a - 0.5 and (b)0.8 for the case with k + Q. The thick solid lines correspond to the growth 
rate of the isobaric mode, Sisobaric- For comparison, the evolutionary rate of the unperturbed central density, apirf), and pressure, 
apirf) are shown by the dashed and dotted lines, respectively. The thin solid lines pointed by arrows correspond to the growth rate 
in the noninteractive mode, Snon-int- 



or <K Lcooi- Since l/k «: Lcoob the sound wave can travel the 
wave length of the perturbation many times during the runaway 
cooling of the unperturbed state. Therefore, the perturbation is 
expected to grow in pressure equilibrium with its surroundings, 
and the mode corresponds to the isobaric mode. An example of 
eigenfunctions (6Q., 6V_t, 6Vy, 6U) for the branch (1) is shown in 
Fig.[3]for rj = 0.75 and k = 8.0. In Fig. [3] it is clearly seen that 
\6n\ <K \6Q.\. This behaviour also implies the isobaric mode. 

The growth rate in the isobaric mode can also be derived an- 
alytically by considering the evolution of a fluid element at the 
centre. The fluid element is assumed to have an isobaric fluctua- 
tion, p - poo(t)+Spoo andP = P()o(0- FromEq. (|4]i, the following 
perturbation equation is obtained: 



dt[-^j^^^-''^'' (r-l)AoP„o Poo 



poo 



From Eqs. (|8]l and (|9]l, we have 

r (r- i)AoPoo ^00 



(1 - W)(fc -f) 
Using Eq. ( |47] ). Eq. ( l46b is rewritten as 



d ( 6pm 
dt\poo 



(2 - a)eoo <Jpo 



(1 - w)(fc - f) Poo 
Equation (l48l l can be integrated to give 

Poo 

Therefore, the growth rate in isobaric mode is given by 



(46) 



(47) 



(48) 



(49) 



1 - 



(2 - a)eoo 

l-O) 

2-a 



r(i - a) 



2-a 



1.72 for rj 
1.20 for 77 = 



r(l - a) I 1.04 for /; = 0.95 



0.10 

0.75 (50) 



In Fig.|2]3, the eigen-value, S, is found to be 1.196 for the case 
with T] - 0.75 and k = 8.0. This is very close to Sisobanc = 1-2 for 
J] = 0.75. From Fig. |2] it is clearly seen that the growth rate ap- 
proaches the corresponding Xisobaric in the large k limit. The ana- 
lytic growth rate is derived by the local argument. Therefore, the 
growth rate is expected to be independent of a global structure 
of the system. iBurkert & LinI (|2000|) performed a linear analy- 
sis on a spatially uniform and isochorically cooling gas. Their 



growth rate in the isobaric mode and our Sisobaiic for 77 = are 
identical, although the spatial structure of the unperturbed state 
is quite different. 

For the density perturbation to grow sufliciently during the 
runaway cooling, it must grow faster than the condensation of 
the cooling layer. This condition can be expressed by Eisobaric > 
ffp. Figure in shows Sisobark, Op and ap as a function of 77 for (a) 
Q'=0.5 and (b) 0.8. From Fig.|4l it is found that Sisobaric is greater 
than Up for all t]. For other a, we can investigate analytically as 
follows: subtracting Op from Sisobaiic, one obtains 



^isobaric ~ 



1 - 
> 1 - ■ 



2-a 

y(l-a) ~ 2~ 
1 



1 



1 + 



2-a 

r(i - a) 



> for ff < 1 . 



(51) 



Therefore, for a < 1, Xisobaiic is more than ap, and the isobaric 
mode can grow in the runaway cooling layer. 



tc - At 



p + Ap 
P - 




tc 




non-interactive 



Fig. 5. Schematic picture of the noninteractive mode. 



(2) The noninteractive mode 

Branch (2) corresponds to the most unstable mode for a- <«; 1 . 
Because the wave length is larger than the cooling length, each 
part evolves independently according to the S-S solutions. We 
call this branch the noninteractive mode. Figure |5] shows the 
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schematic picture of the noninteractive mode. In Fig. |5j simi- 
larity variables have an initial fluctuation. For example, we con- 
sider two different regions, "A" and "B", where pA - p + and 
Pb - p. Due to the difference of Ap, the regions "A" and "B" 
have different collapse time, t^-At and t^, respectively. Omitting 
any terms that do not grow, we find the time evolution of differ- 
ence, Ap, to be 



Ap 
Po(0 



= -Q(^)Af 



1 dlnQ 



On + 



1 - w dln^ 



1 



tr - t 



(52) 



Therefore, in the zooming coordinate, the density perturbation 
grows as ^Q(^) (t^. - Other perturbed variables also grow 
in the same power law. Therefore, comparing with Eq. ( |27] i, the 
growth rate is given by 



1, 



(53) 



which is independent of a and t]. The noninteractive mode arises 
from the fluctuation of the collapse time, t^, due to density and 
pressure perturbations. From the physical mechanism, the per- 
turbation grows in the same way as the unperturbed state. In 
other words, the perturbation of the isobaric (isochoric) cooling 
layer grows isobarically (isochorically). 

We investigate whether the noninteractive mode grows suf- 
ficiently during the runaway cooling or not. First, we consider 
the case with rj > T/eq- Since the perturbation grows isobarically, 
2^non-int IS Compared to ap = r]/(2 - a). From Fig. 21 it is seen 
that the growth rate, Enon-int, is higher than Op for all t] > rj^q - 
0.75 and 0.86 for a - 0.5 and 0.8, respectively. Therefore, the 
noninteractive mode can grow sufficiently. Next, we consider the 
cooling layer with 77 < //eq. In this layer, since the perturbation 
grows isochorically, Xnon-int is compared to = (1 - ri)/{l - a). 
Analytically, it is found that the pressure perturbation can only 
gro w for T] > a. 

iKovama & Inutsukal (12000') performed a linear analysis of a 
spatially uniform and isobarically cooling gas in their appendix. 
However, they did not find the noninteractive mode. This is 
because they fixed the collapse time to be spatially constant, 
and it was a s sumed not to be influenced by perturbation. 
iBurkert & LinI (l200d) performed linear analysis of a spatially 
uniform and isochorically cooling gas by taking account of the 
time evolution of the unperturbed state. They showed that a 
perturbation cannot grow in the condition for long wave length 
limit. Our result is consistent with theirs. 

(3) The shear mode 

For /(■ <s; 0, there is a solution where physical variables 
are very small except for 6Vy which is spatially constant, and 
the eigenvalue is cr - u. The similar mode was found by 
iMcNamara ' ('1993|, w ho inv estigated Tl of a uniform granular 
medium. McNa mara] (1 19931) called this mode the shear mode. 
The growth rate is given by 




for 77 = 0.10 
for 77 = 0.75 
for 77 = 0.95 



(54) 



In Fig. 12] it is seen that each growth rate in branch (3) has the 
corresponding value of Sshear for k <^ \. The physical meaning 
of this mode can be understood as follows: for a: <K 1, since the 
effect of the pressure gradient with respect to y is very weak, 
the gas can freely stream with almost constant velocity, y,,, in 
the y direction. On the other hand, the central sound speed, 
coo(f)> decreases as oc (t^. - f)"^^'"'^*. Therefore, the ratio of the 
dynamical velocity to the thermal velocity, Vy/coo(f)^ grows with 



time as oc {t^ - f)""^('"'^*, indicating that the growth rate is given 
by Eq. ( |54b . For the case with larger wave number, the effect of 
the pressure gradient becomes important. Therefore, the fluid 
element cannot stream freely, and the growth rate is lower as 
shown in Fig.|2] 

(4) The free-streaming mode 

For large k, there is another mode in which the velocity per- 
turbation in the x-direction is much greater than that in the y- 
direction, \6Vxi)\ » I^Vyol. We call this mode the free-streaming 
mode. From Eq. ( |45] | with 5Qo = OTo = Slim - 0, we obtain 



■-free 



2yg 

1 -O) 




for 77 = 0.10 
for 77 = 0.75 
for 77 = 0.95 



(55) 



In the free-streaming mode, the growth of the velocity perturba- 
tion in the .jc-direction is hampered by the pressure gradient of 
the unperturbed state. Therefore, the growth rate is less than the 
shear mode. 

(5) ^ = mode 

The growth rate in this branch for k coincides with the case 
with k -Q, which is obtained in Sect. [3] 



4.2. Linear analysis considering tine time evoiution of k 

The static approximation is valid only if S is much larger 
than |dln/c/dlnft|. However, from Fig. |2l it is found that the 
growth rate of the most unstable mode for each k is smaller 
than |dln/c/dlnft| whose values are 1.93, 1.5, and 1.37 for 
77 = 0.1, 0.75, and 0.95 with a = 0.5, respectively [see Eq. ([39]ll. 
Therefore, in this section, we perform a Unear analysis consider- 
ing the time evolution of k{t) using direct numerical integration. 
The upwind difference method is used as the numerical method 
to solve the perturbation equation. The region of calculation is 
from ^ = to ^ = 100 in the zooming coordinate. 

We impose the boundary conditions at f = and ^ = 100. 
The even mode is set as the boundary condition at ^ = 0. The free 
boundary condition is set at ^ = 100, but it does not influence 
the inner region since the gas flows out supersonically from the 
outer boundary of the zooming coordinate. As an initial state, we 
adopt the eigenfunction of the isobaric mode for k - 30 which is 
obtained in Sect. 14.11 

By solving Eqs. (l2TI)-(l24b. a time evolution of dQ is obtained. 
During the calculation, the growth rate at r is evaluated by 



1 



1 - (jj dr 



{ln(5Q(^ = 0,T)) 



(56) 



at each instant of time. For comparison with the result of the 
static approximation, we focus on a relation between /c(f) and the 
growth rate of the density perturbation at the centre, Snum. Figure 
|6]shows the growth rate as a function of k{t) at each instant 
of time for a - 0.5 with (a)?7=0.1 and (b)0.75. The nondimen- 
sional wave number, k(t), decreases with time as shown in Eq. 
(l37T i. Therefore, in Fig.|6] the direction of time is from the right 
to the left. For comparison, the approximate dispersion relations 
of branches (l)-(2) in Fig. [2] are superimposed by the dashed 
lines. In both of Figs. |6^ and|6]3, the behaviour of the growth rate, 
^^nuiTi, moderately agrees with that of the approximated disper- 
sion relations. For tc » 1, or initial phase, the growth rate agrees 
with Sisobaric- This is because the growth rate does not depend 
on K in the short wave length limit. As k decreases and reaches 
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K K 

Fig. 6. Growth rate obtained by results of numerical linear analysis for rj- 0.1(a) and 0.75(b). The thick and dashed lines indicate 
the results of numerical linear analysis and those of approximate linear analysis. 



about 1, Snum begins to decrease. For k ^ I, l,„um approaches 
asymptotically Snon-im, where the effect of k is negligible. The 
effect of time-depending k is notable only for 0.1 < a: < 10. 
Smoother dependence of the growth rate on k is obtained than 
the approximate dispersion relation. 

5. Discussion 

5.1. Astrophysical implication 

In this section, we estimate the fragmentation time of the cool- 
ing layer formed by a collision of WNM. We consider a head- 
on collision between two thermally equilibrium gases in which 
density and pressure are Pwnm - O.Slmn cnT^ and Pwnm - 
3.5 X 10^A:b dyne cm"^. Two clouds collide along the x-axis at 
f = and X - with velocity 20 km s i.e., the mach number 
is 2.17. The c ooling function in this temperature region fitted by 
iKovama & In utsuka (2002) as follows: 



10= 



pX. - «(-r + nA) erg cm s 
where 

r = 2 X 10"^^ 



A 
f 



l.OxlO'^ exp 



118400 \ , , ,^ 
+1.4x10 

r + 1000/ 



-2 



exp 



92 
'Y 



(57) 



(58) 



,(59) 



where n is the number density of the gas. We perform numeri- 
cal calc ulation using the one-dimensional Lagrangian Godunov 
scheme dvan Leeiill997 ') between x - and x - - 3.26pc. 
At x = and x — Lj^, we adopt the reflecting and free bound- 
ary conditions, respectively. As the initial condition, we add the 
following density fluctuation: 



pit -Q,x) - 



Pwnm 



»^max \ '-'X 



(60) 



where we set /max = 10, A = 0.5, and the phase, 0,, is given by 
random number between and 2;7r. The maximum and minimum 
number density at the initial state are 0.72 cm"-' and 0.44 cm"-', 
respectively. 

Figure |2] shows the evolutionary track of the gas at the cen- 
tre, X - Q. Due to the shock compression, the temperature of 
gas increases suddenly, and the postshock region enters the ther- 
mally unstable phase from the initial stable state. The TI leads 



CO 




cm 




^ 10" 


o 


II 

><_ 




Q. 





10" 



Runaway Cooling 



10 




n(x=0) [cm"3] 



Fig. 7. Evolutionary track of the gas at the centre after the head- 
on collision (the thick solid line). The ordinate and abscissa axes 
indicate the pressure and the number density, respectively. The 
thin solid line indicates the thermally equilibrium curve. The 
filled circle represents the gas at the preshock region. 




0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 
t [Myr] 

Fig. 8. Time evolution of anum (the solid line) and the cool- 
ing length (the dotted line), which are evaluated at the centre . 
The dashed line indicates the time evolution of the critical wave 
length defined in Sect. | 



to the cooling layer condensing in a runaway fashion. The cool- 
ing length in Fig. |8] rapidly decreases with time until it reaches 
minimum value 1.5 x 10"^ pc at f = 0.917 Myr. During runaway 
cooling, the gas condenses isobarically until it reaches the stable 
CNM phase, although some pressure oscillation is seen in Fig. 
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x[pc] xipcKi-t/y-""*''"' 

Fig. 9. Time evolution of (a)the number density and (b)the temperature, respectively. The rescaled (c) number density, nfj^'^'"^ and 
(d) temperature, 7'f('"°'*('"''^/'''(2""'^l as a function of the rescaled coordinate, where fc and rj are set to 0.915 Myr and 0.98, 

respectively. The thick lines in (c) and (d) indicate the corresponding S-S solution. 



|2] We focus on the property of this runaway condensing layer 
First we evaluate Q'niim(0^ which is defined by 

Q'num = ;3T, (61) 

at each instant of time. Figure [8] shows the time evolution of 
Q^num- During 0.4 < f/Myr < 0.9, it is seen that anum ~ 0.61 
is approximately constant. Therefore, in this period, the flow is 
expected to be approximated well by the S-S solutions. 

Figures |9] show the time evolutions of the distributions of 
(a) the number density and (b) the temperature at f = 0.72, 
0.80, 0.88, and 0.90 Myr, respectively. The rescaled number den- 
sity «?, and temperature, yf, are shown m 

Figs. |9}; and|9}J as a function of rescaled coordinate, xf» 
where fc, rj, and a are set to 0.915 Myr, 0.98, and 0.61, respec- 
tively. From Figs.|9}; and|9}i, it is clearly seen that the S-S solu- 
tion (a - 0.61,7; - 0.98) describes the results of the numerical 
calculation well. There are two reasons they are described by 
the S-S solution well. One is that Q'num is almost constant, and 
the cooling function is approximated well by Eq. ([1]) during the 
runaway condensation. The other reason relates to the stability of 
S-S solutions. The one-dimensional calculation solves the evo- 
lution only in the direction parallel to the condensation. From 
the linear analysis in Sect. [3] such a perturbation with k - Q can- 
not grow enough during the runaway condensation. That is why 
the S-S solution is realized even though the density is initially 
fluctuated. 

However, the result is expected to be qualitatively different 
if the perturbation with k i^Q exists. The perturbation with k i^Q 
on the isobarically runaway condensing layer grows as oc t^^ (see 
Sect.Hl). The growth rate is independent of wave numbers. From 
Figs.[8]and|9l the gas evolves obeying the S-S solution between 



t - 0.4 and 0.9 Myr. If the layer has a perturbation ik + 0) with 
amplitude dpo at f = 0.4 Myr, the perturbation grows as 

&p{t) ^ 1 - 0.4Myr/fe ^g^) 
(5p() 1 - f/fc 

The epochs when the perturbation grows by a factor of 10 and 
100 are t - 0.86 and 0.91 Myr, respectively. At 0.86 Myr, the 
central number density is still as low as 38 cm"-'. Therefore, 
the condensing layer is expected to fragment quickly, and CNM 
clu mps will form. 

llnoue & Inutsukal (l2008l) investigated TI in a shock com- 
pressed region formed by WNM-WNM collision using two- 
dimensional, two-fluid magnetohydrodynamical simulation. The 
initial condition there is the same as that in our one-dimensional 
simulation except for the dimension and the way to add density 
fluctuation in the WNM. Our linear analysis is applicable to the 
gas evolution before CNM clouds form in their unmagnetized 
case. In the initial stage of their simulation, in the postshock re- 
gion, thermally unstable gas initially condenses in a layer struc- 
ture. Around t - 0.5 Myr when the region of the highest density 
reaches CNM stable state, both sm all CNM clou dlets and fila- 
mentary CNM clouds are generated dlnou According to 
our linear analysis, an isobarically condensing gas layer is ex- 
pected to split into fragments that have a variety of lengths in the 
transverse direction. Therefore, our linear analysis agrees with 
their simulation qualitatively. 

After the first CNM clouds formation, the mass of these 
CNM clouds continues to increase by the accretion of unsta- 
ble gas by the shocks. Moreover, some CNM clouds coalesce 
into larger clouds and o t hers f ragment into smaller ones by 
the turbulent flow dlnoud l2009l) . Therefore, the size and mass 
of CNM clouds varies with time. To understand this phase, 
which is beyond the scope of this paper, a statistical approach 
is probably needed. A statistical theory has been proposed by 
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iHennebelle & Audit! (l2007h using Press-Schechter formalism 
dPress & SchechteR 1 19741) . assuming that the CNM clumps are 



generated from density fluctuations within WNM whose power 

spectrum is assumed to be Kolmogorov. 

Comparing timescales. lHeitsch. Hartmann & Burkert! (12008!) 
discuss the effect of TI, the nonUnear thin-shell instability, 
Kelvin-Helmholtz instability, and gravitational instability in var- 
ious densities, temperatures, and scales of fluctuations. In their 
paper, it is assumed that the gas is unstable only if the scale of 
the fluctuation is smaller than the sound crossing scale. However, 
from our linear analysis, perturbations can grow regardless of 
their scales as long as they are flat rather than spherical. 

5.2. The growth rate fori < a <2 

Although the linear analysis on the S-S cooling layer is limited 
for or < 1, the thermal stability of the gas for o- > 1 is also 
roughly understood from Balbus's criterion. For 1 < a < 2, 
the gas is isobarically unstable, but it is isochorically stable. For 
a > 2, the gas is thermally stable. In this section, we investigate 
the stability of the gas for 1 < a < 2 during cooling within the 
large and small scale limits. 

5.2.1. The isobaric mode 

For the case with small wave length, perturbation is expected 
to grow isobarically. By comparison of our results with previ- 
ous studies in the literature, it is found that the growth rate in 
the isobaric mode is independent of the global structure of the 
unperturbed state. Therefore, from local arguments, the growth 
rate in the isobaric mode of the gas with 1 < a < 2 can also be 
estimated. 

As an unperturbed state, we adopt a cooling gas element 
whose scale is assumed to be much smaller than the cooling 
length. In this case, the element cools isobarically. From Eq. (0]), 
the time evolution of the unperturbed gas is given by 



Comparing Eq. d68T l with Eq. d63T l, one can see that the perturba- 
tion grows more slowly than the unperturbed state for I < a <2. 
Therefore, the gas is expected to be difficult to fragment during 
runaway cooUng if 1 < a <2. 

5.2.2. The noninteractive mode 

A cooling layer that evolves isobarically is considered. The time 
evolution of the central density is the same as Eq. ( l63T l. When 
the scale of perturbation perpendicular to the condensation is 
too large to interact with other regions, each region evolves in- 
dependently. Here, we focus on the time evolution of density 
perturbation at the centre {x - 0). Initial fluctuation of the cen- 
tral density, dp,, creates the fluctuation of the cooling time, Af. 
The relative amplitude of density perturbation at x = is given 
by 



6p 1 

-^--(p,+5p,) 1-- 
po Po \ t__, - At 



cool 



-I/(2-ff) 



- 1. 



(69) 



Linearizing Eq. 
have 

6p _ 1 Af 
PO 2- a t' 



with omitting terms that do not grow, we 



(70) 



cool 



-U(2-a) 



1 



pit) = p,. 1 1 - — I , — = (2-a)y"-Uj-l)Pr'pJ-", (63) 



cool 



where p, and P, represent the initial density and pressure, respec- 
tively. In the above unperturbed state, we consider the following 
isobaric perturbation: 



Comparing Eq. ( |70] i with Eq. (|63j, we can see that the perturba- 
tion grows more slowly than the unperturbed state for 1 < a < 2. 
Therefore, the gas is expected to be diflicult to fragment for 
1 < < 2 for the large-scale perturbation, as well as the small 
scale. 

5.3. Effects of thermal conduction 

In this paper, the thermal conduction is neglected for simplic- 
ity. However, for large wave number, the therm al conductio n 
is expected to stabilize Tl in the cooling layer dFieldl 11965!) . 
Therefore, there is a critical wave number, kcm, such that per- 
turbation with a larger wave number is stabilized by the thermal 
conduction. 

First, we evaluate k^^i-it using an order estimation. Using the 
characteristic time scale of the thermal conduction, fdiff , the dif- 
fusion equation is given by 



p = p()(f) + Spit), 
and 



(64) _L 



00 



fdiff \r- 1 



k^K{TQo)T^ 



00, 



(71) 



(65) 



where K is the thermal conduction coefficient. From Eq. ( fTTT i. 
the diffusion timescale is given by 



where subscript "0" indicates the unperturbed state, and 6p is the 
density perturbation. Linearizing Eq. (|4|i, one obtains 



fdiff - 



00 



r,-2 



(r- l)KiToo)Ti 



(72) 



00 



d 

d^' 



(2- 



6p 
[po 

Using Eq. (|63]), Eq 
1 



a)yr-Vo-"Po-'-- 
Po 



(66) 



d 

d^' 



IS rewritten as 

I 



From Eq. d72] i. one can see that the diffusion timescale is small 
for large wave number If fdiff < fcoob the thermal conduction is 
expected to stabilize Tl. Therefore, ka-it can be derived on the 
condition fdiiF ~ fcooi as 



\P0l 'cooll 'cool/ 

Equation ( |67] l is easily integrated to give 

-1 



(67) k,,n = 



-'cool 



where = 



Poocoo 
(y-l)KT, 



(73) 



00 



Sp 
Po 



— oc 1 



t 

'cool 



(68) 



iField! (Il965h derived similar critical wave number to Eq. ( l73T l. 
Since the unperturbed state is time dependent, A;ciit also evolves 
with time. Detailed evolution of kan depends on K.lnT < 6000, 
we adopt K ^2.5x10^ yff ergs cm ' K^' s^' (iParkerL[T953l) . In 



Kazunari Iwasaki and Tom Tsuribe: Fragmentation of a dynamically condensing radiative layer 



11 



this case, from Eq. (fTTt . the time evolution of k^rit can be derived 
analytically as 

dln/^crit _ {2a - 1)?? - (2 - g)(3 - 2a) 
din ~ M2-a)(l - a) 

< for 7/ = 1 

< for 7/ = 



4(2 - a) 
3 -2a 

'4(1 - a) 



(74) 



For 77 = 1 and 77 = 0, it is found that Eq. ( |74] | is negative for a < 
1. Since Eq. (l74l i is the linear function for 77, dlnfccrit/dlnf, is 
negative for all rj. Therefore, kah increases with time. This means 
that an initially stable perturbation with wave number (k > kcnt) 
becomes unstable at a certain epoch when k^-iyi catches up with k. 

Figure|8]shows the time evolution of the critical wave length. 
Adit - 2jT/kciit- In Sect. 15. II the time evolution of the condensing 
layer can be described by the S-S solution with (anum = 0.61, 
77 = 0.98). In Fig. [8] we see that Aa-it decreases with time during 
the runaway condensation. Figure|8]also shows that Acdt < icooi, 
or Kait = kcntLcoo\ > 1- This means that the effect of thermal 
conduction on the dispersion relation (Fig.|6]l always appears in 
the isobaric regime during the runaway condensation. 

6. Summary 

In this paper, we have investigated the stability of S-S solutions 
describing the runaway cooling of a radiative gas by linear anal- 
ysis. The results of our investigation are summarized as follows, 

1 . For the case with perturbation only parallel to the flow (k - 
0), the S-S solutions are unstable. However, the growth rate 
is too low to become nonlinear during the runaway cooling. 
Actually, the S-S solutions are realized in one-dimensional 
hydrodynamical calculations in IT08 and in Sect. l5.ll in this 
paper. 

2. For the case with transverse perturbation (k + 0), there are 
several unstable modes in the S-S solutions. The most unsta- 
ble modes are the isobaric mode for » 1 and the noninter- 
active mode for k-^ 1 . In the isobaric mode, the perturbation 
grows in pressure equilibrium with its suiToundings. On the 
other hand, the noninteractive mode is originated from each 
region in the layer condensing independently. Under a static 
approximation, we derive the approximated dispersion rela- 
tion. The results of direct numerical integration of the time 
evolution agree with those using the static approximation. 

3. The S-S solutions for 77 > o' are unstable for any wavelength. 
Especially, if the unperturbed state is isobaric, the growth 
rate is independent of wave number. Therefore, fluctuations 
in various scales grow simultaneously, and the gas layer is 
expected to split into fragments with various scales. The S-S 
solutions for 77 < c are only unstable in the isobaric mode. 

Our linear analysis predicts that the cooling layer splits into 
fragments quickly even if the size is greater than the local cool- 
ing length. Our linear analysis is qualitatively consistent with the 
results of recent multi-dimensional numerical simulations until 
CNM clouds form, but the evolution of CNM clouds is beyond 
the scope of this paper. 
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Appendix A: Supplements for derivation of basic 
equations in the zooming coordinate 

Here, we present preparation for derivation of basic equations 
(fT0li-(fT2li in the zooming coordinate given by Eq. (|5]l. In the 
zooming coordinate, the physical variables, Q(t,x,y), are given 
by the following unified form: 



ql{\-u>) 



(A.l) 



where Q{t,x,y) corresponds to [p, Vi, Vv,P] (see Eq. |6]l, and 
= [O, Vx, Vy, n] are the physical variables in the zooming co- 
ordinate. From Eqs. Q-®, a parameter, q, is given by 



(-P for Q^p 
q-<-(ji) fox Q- Vx, Vv . 
-2a) -p for Q^P 



(A.2) 



The temporal and spatial derivatives of Q{t, x, y) in the ordinary 
coordinate can be expressed in the zooming coordinate as 



dQ\ 
dt I 



dQ 

dx 



^e,e„(o__.eo(o(|)^- 



vo(0 



Go(f) 



xoit) 
xo(f) 



aw?, 

dy 



(A.3) 
(A.4) 
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respectively, where we use 
vo(0 xo(t) 1 



(A 5) 

xo(t) xo(t) (l-w)fcf. 

from Eq. (|7]i. Using Eqs. ( IA.3l l and ( IA.4l i. the Lagrangian time 
derivative of Q is given by 



(d_ 



d 



d 



+ vo(f)V,— + vo(f)Vv— e(f, x,y) 



dx 



dy 



vo(t) 



Qo(t) 



d d d 

q + ^ + (v, + ^) — + v,xo(t)— 

OT of oy 



0(T,f,y).(A.6) 



xo(0 

Using Eqs. (|A3]i-(|a;6]i. 

^=r-'AoMO^-"^o(^r', and P,it) ^ P^^^AA.l) 
xo(t) y 

basic equations (fT0]i- (fT2l i are derived in the zooming coordinate. 



Appendix B: Detailed expression of 

In this appendix, we provide the detail expression of A,^ as 



An =- 



Vo j 



Xl 



Xl 



^_^(lnno)' + --!i(a-2)eo 

r Vo 

An = -Vo(lnQo)' + cr + w + - _2-(lnnon7)'. 



^0 

An = — 



(In no) 



, o- + (a-l)yei) 



Vo 



Ai4 = -^.xo(f)Vo, 



(B.l) 

(B.2) 

(B.3) 
(B.4) 



A21 = ^(Vo(lnno)'-(ff-2)reo), 

r 



(B.5) 



A22 = X2(ln Qo)' - Vo(tr + oj + V'q) + -^(In HoQ^)', (B.6) 

An = -V0A14, (B.7) 
A24 = ^^0(0^0' (B-8) 



A31 = -X2(lnno)' + (ff - 2)yoreo, 

A32 = r {- Vo(ln Qo)' + cr + + - -^(In OoQ^) 

A33 = X^anDo)' - Vo{cr + (a - l)reo), 
A34 = -y^xo(OVo, 



(B.9) 

(B.IO) 

(B.ll) 
(B.12) 



A41 = 0, 
A42 = 0, 

A43 ■■ 
and 



- x-^ 
^0 



Vo 



(w + cr), 



X^ 

AAA^kxQ{t){Vl-Xl)^. 

7 



(B.13) 
(B.14) 

(B.15) 



(B.16) 



